#include "qr.h"
#include "householder.h"
#include <cmath>

#include <iostream>
using namespace std;

QR::QR(Matrix &matrix, double error): error(error), matrix(matrix){
	this->x = Matrix(matrix.getRow(),matrix.getColumn(),IDENTITY_MATRIX);
	this->lambda = Matrix(matrix.getRow(),matrix.getColumn(),IDENTITY_MATRIX);
}

QR::~QR(){
}

void QR::calculate(){
	
// 	std::cout << "matrix original" << std::endl;
// 	matrix.print2();
	
	x = Matrix(matrix.getRow(),matrix.getColumn(),IDENTITY_MATRIX);	//  x - auto vetores
	
	//utilizando householder
	Householder method;
	matrix = method.calculate(matrix, x);
	std::cout << "feito Householder" << std::endl;
	
	
	//metodo de qr em si
	Matrix current = Matrix(matrix);
	Matrix r = Matrix(matrix);	//  r - triangular superior
	Matrix jt, qt;			//  produtorio(jt[i][j]) = qt
	double theta = 0.0, acum = 0.0;
	//double acum_;
	
	do{//iteracao k
		//iniciando sempre a Qtransposta do passo k com a identidade
		qt = Matrix(matrix.getRow(),matrix.getColumn(),IDENTITY_MATRIX);
		
		for (unsigned int j = 0; j < matrix.getColumn() - 1 ; j++){		// coluna
//  			for (unsigned int i = j + 1; i < matrix.getRow(); i++){ 	// linha
				unsigned int i = j + 1;
				//iniciando smepre a matriz de jacobi com a identidade
				jt = Matrix(matrix.getRow(),matrix.getColumn(),IDENTITY_MATRIX);
				
				//calculando teta
				//theta = acos( (current[j][j])/sqrt( pow(current[i][j],2) + pow(current[j][j],2) ) );
				theta = std::atan(current[i][j]/current[j][j]);
				
				//preenchendo a matriz de jacobi como matriz de rotação
				jt[i][i] = cos(theta);
				jt[j][j] = cos(theta);
				jt[i][j] = -sin(theta);
				jt[j][i] = sin(theta);
				
				//calculo do r
				r = jt * r;
				
				//calculo do q transposto
				qt = jt * qt;
//  			}
		}//ao final tem o rk e o qtk
		
		//Ak+1 = rk*qk
		current = r * qt.transpose();
		
		//preparando r
		r = current;
		
		//atualizando a matriz dos autovetores
		x = x * qt.transpose();
 	//	acum_ = acum;
		acum = 0.0;
		for (unsigned int j = 0; j < matrix.getColumn() - 1 ; j++)	// coluna
			//for (unsigned int i = j+1; i < matrix.getRow(); i++)	// linha
					acum+= current[j+1][j];
		lambda = current;
//		cout << acum << endl;
	}while(std::abs(acum)>this->error/* && std::abs(acum_ - acum) > this->error*/);
}

void QR::printLambda(){
	lambda.print();

}

void QR::printVector(){
	x.print();
}

Matrix QR::getLambda(){
	return lambda;
}

Matrix QR::getVector(){
	return x;
}
